## This code manages data and conducts the analysis
## reported in Ahlquist, John S. 2010. "Building and Using Strategic Capacity"
## American Political Science Review
## Note that S+ code for WSEV is available from http://polmeth.wustl.edu/polanalysis/vol/10/heagerty.window-var.ssc
##
## Last updated 021510
######################################

library(geepack)
library(xtable)
library(MASS)
library(lattice)
library(zoo)

#### Functions used later in the analysis ####

## Heagerty et al. WSEV code, slightly modified to work in R

window.var<-function(y, X, fit, m, tvec, step = 1){
#
# PURPOSE:  compute window resampling variance estimates for
#           binary spatial time series
#
#
##
###  Let "m" be length of the window
##
#
#  m <- trunc( C * N^{d/(d+2)} )
#
#
##
###  Let N be the number of time points
##
#
    N <- length(unique(tvec))
    nobs <- length(y)   #
##
### compute number of sub-windows to consider
##
#
    n.windows <- trunc((N - m + 1)/step)    #
##
### preliminary calculations
##
#
    antilogit <- function(x)
    {
        exp(x)/(1 + exp(x))
    }
#
    p.fitted <- fitted(fit)
    v.fitted <- p.fitted * (1 - p.fitted)
    R <- (y - p.fitted)
    U.N <- t(X) %*% R
    H.N <- t(X) %*% (as.vector(v.fitted) * X)    #
##
### compute U U^T and the one-step U U^T 
##
# 
    p <- ncol(X)
    eUU <- matrix(0, p, p)
    eUU.1 <- matrix(0, p, p)
    mU <- matrix(0, p, 1)   #
    utimes <- unique(tvec)
    utimes <- utimes[order(utimes)]
    for(j in 1:n.windows) {
        start <- utimes[(j - 1) * step + 1]
        stop <- start + m - 1   #
        keep <- (tvec >= start) & (tvec <= stop)    #
        X.j <- X[keep,  ]
        U.j <- t(X.j) %*% R[keep]   #
        H.j <- t(X.j) %*% ((v.fitted[keep]) * X.j)
        beta.1 <- fit$coef + as.vector(solve(H.N - H.j) %*% (U.N - U.j))    #
        lp.j <- as.vector(X.j %*% beta.1)
        U.j.1 <- t(X.j) %*% (y[keep] - antilogit(lp.j)) #
        eUU <- eUU + ((N/m) * (U.j %*% t(U.j)))/n.windows
        eUU.1 <- eUU.1 + ((N/m) * (U.j.1 %*% t(U.j.1)))/n.windows
        mU <- mU + U.j.1/n.windows
    }
#
    I.N <- solve(H.N)   #
    center <- T
    if(center)
        eUU.1 <- eUU.1 - (N/m) * mU %*% t(mU)
    var.est.0 <- I.N %*% eUU %*% I.N
    var.est.1 <- I.N %*% eUU.1 %*% I.N  #
    out <- list(title = "Window resampling for estimating functions", 
        m = m, step = step, n.windows = 
        n.windows, independence = I.N, 
        window.0 = var.est.0, window.1 = var.est.1) #
    out
}
## End WSEV code


## A hacked version of apsrtable ##
table.edit<-function (..., se = c("robust", "vcov", "both"), model.names = NULL, 
    model.counter = 1, digits = 2, stars = 1, lev = 0.05, align = c("left", 
        "center", "right"), order = c("lr", "rl", "longest"), 
    omitcoef = NULL, Sweave = FALSE, Minionfig = FALSE) 
{
    x <- character(0)
    signif.stars <- TRUE
    order <- match.arg(order, c("lr", "rl", "longest"))
    opts <- match.call(expand.dots = FALSE)
    se <- match.arg(se, c("robust", "vcov", "both"))
    align <- substr(align, 1, 1)
    se = match.arg(se, c("robust", "vcov", "both"))
    align <- match.arg(align, c("l", "c", "r"))
    adigits <- ifelse(align == "c", -1, digits)
    models <- list(...)
    nmodels <- length(models)
    if (!Sweave) {
        x <- cat("\\begin{table}[!ht]\n\\caption{}\n\\label{}\n")
    }
    if (Minionfig) {
        x <- cat("%Uncomment the following line and the end one to change figure versions\n%if you are using a full-featured family such as Minion Pro.\n\\figureversion{tabular}\n")
    }
    model.summaries <- lapply(models, function(x) {
        s <- summary(x)
        if (!is.null(x$se) && se != "vcov") {
            est <- coef(x)
            if (class(x$se) == "matrix") {
                x$se <- sqrt(diag(x$se))
            }
            s$coefficients[, 3] <- tval <- est/x$se
            s$coefficients[, 4] <- 2 * pt(abs(tval), length(x$residuals) - 
                x$rank, lower.tail = FALSE)
            s$se <- x$se
        }
        return(s)
    })
    if (is.null(model.names)) {
        m.first = model.counter
        m.last = m.first + (nmodels - 1)
        model.names = paste("Model", m.first:m.last)
    }
    else if (!is.null(model.names) && (length(model.names) < 
        nmodels)) {
        m.first = length(model.names) + 1
        model.names = c(model.names, paste("Model", m.first:nmodels))
    }
    mlength <- sapply(model.summaries, function(x) length(coef(x)))
    longest <- which.max(mlength)
    if (order == "rl") {
        modelorder <- nmodels:1
    }
    else {
        modelorder <- 1:nmodels
    }
    if (order == "longest") {
        coefnames <- names(coef(models[[longest]]))
    }
    else {
        coefnames <- names(coef(models[[modelorder[1]]]))
    }
    for (i in modelorder) {
        matched <- match(names(models[[i]]$coef), coefnames, 
            nomatch = 0)
        unmatched <- which(is.na(matched) | matched == 0)
        coefnames <- c(coefnames, names(models[[i]]$coef)[unmatched])
    }
    incl <- rep(TRUE, length(coefnames))
    names(incl) <- coefnames
    if (!is.null(omitcoef)) {
        incl[omitcoef] <- FALSE
    }
    model.summaries <- lapply(model.summaries, function(x) {
        pos <- match(rownames(coef(x)), coefnames)
        x$var.pos <- pos
        return(x)
    })
    out.table <- lapply(model.summaries, function(x) {
        var.pos <- x$var.pos
        model.out <- model.se.out <- star.out <- rep(NA, length(coefnames))
        model.out[var.pos] <- x$coefficients[, 1]
        #star.out[var.pos] <- apsrStars(x$coefficients, stars = stars, 
        #    lev = lev, signif.stars = TRUE)
        model.out <- ifelse(!is.na(model.out), paste(formatC(model.out, 
            digits = digits, format = "f")), "")
        model.se.out[var.pos] <- x$coefficients[, 2]
        if (!is.null(x$se) & se %in% c("robust", "both")) {
            model.se.out[var.pos] <- x$se
        }
        model.se.out <- ifelse(!is.na(model.se.out), paste("(", 
            formatC(model.se.out, digits = digits, format = "f"), 
            ")", sep = ""), "")
        if (se == "both" && !is.null(x$se)) {
            model.se.out[var.pos] <- ifelse(model.se.out != "", 
                paste(model.se.out, " [", formatC(x$coefficients[, 
                  2], digits = digits, format = "f"), "]", sep = ""), 
                "")
        }
        model.out <- rep(model.out[incl], each = 2)
        model.se.out <- rep(model.se.out[incl], each = 2)
        pos.se <- (1:length(model.out))[(1:length(model.out)%%2 == 
            0)]
        model.out[pos.se] <- model.se.out[pos.se]
        model.info <- list("$N$" = formatC(sum(x$df[1:2]), format = "d"), 
            "$R^2$" = formatC(x$r.squared, format = "f", digits = digits), 
            "adj. $R^2$" = formatC(x$adj.r.squared, format = "f", 
                digits = digits), AIC = formatC(x$aic, format = "f", 
                digits = digits), BIC = formatC(((x$aic - 2 * 
                (length(x$coef))) + log(sum(x$df[1:2])) * length(coef(x))), 
                format = "f", digits = digits), "$\log L$" = formatC(((x$aic - 
                2 * (length(x$coef)))/-2), format = "f", digits = digits))
        attr(model.out, "model.info") <- unlist(model.info)
        return(model.out)
    })
    out.matrix <- matrix(unlist(out.table), length(coefnames[incl]) * 
        2, nmodels)
    out.matrix <- cbind(rep(coefnames[incl], each = 2), out.matrix)
    out.matrix[(row(out.matrix)[, 1]%%2 == 0), 1] <- ""
    #out.info <- matrix(unlist(lapply(out.table, attr, "model.info")), 
    #    6, nmodels, dimnames = list(names(attr(out.table[[1]], 
    #        "model.info")), NULL))
    #out.info <- as.matrix(out.info[apply(out.info, 1, function(x) sum(x == 
    #    "") != nmodels), ])
    #out.info <- cbind(as.character(rownames(out.info)), out.info)
    #out.matrix <- rbind(c("%", model.names), out.matrix)
    #out.matrix <- rbind(out.matrix, out.info)
    #out.matrix[, -1] <- format(out.matrix[, -1])
    #out.matrix[, 1] <- format(out.matrix)[, 1]
    out.matrix <- rbind(c("", paste("\\multicolumn{1}{", align, 
        "}{", model.names, "}", sep = "")), out.matrix)
    out.matrix <- apply(out.matrix, 1, paste, collapse = " & ")
    x <- cat(paste("\\begin{tabular}{", align, paste("D{.}{.}{", 
        rep(adigits, nmodels), "}", sep = "", collapse = ""), 
        "}", sep = ""))
    x <- cat("\\hline\n")
    x <- cat(paste(out.matrix, collapse = "\\\\ \n"))
    se <- ifelse((se != "vcov" && sum(unlist(lapply(models, function(x) !is.null(x$se))) > 
        0)), "robust", "vcov")
    cat("\\\\\\hline\n\\multicolumn{", nmodels + 1, "}{l}{", 
        ifelse(se != "vcov", "Robust s", "S"), "tandard errors in parentheses}", 
        ifelse(se == "both", paste("\\\\\n\\multicolumn{", nmodels + 
            1, "}{l}{", "Na\\\"ive standard errors in brackets}", 
            collapse = "", sep = ""), "\\\\"), sep = "")
    cat("\n\\end{tabular}\n")
    if (Minionfig) {
        cat("\n\\figureversion{proportional}\n")
    }
    if (!Sweave) {
        cat("\n\\end{table}\n")
    }
    return(invisible(x))
}

## END table.edit

### END helper functions


#### Analysis and graphics ####
#dowload the file AhlquistAPSR2010data.csv and set R's working directory accordingly

data<-read.csv("AhlquistAPSR2010data.csv")

data$enpp<-log(data$enpp)
data$pop<-log(data$pop)

#setting up dummies for coordination
Level.gt1<-Level.gt2<-Level.gt3<-Level.gt4<-Wcoord.gt1<-Wcoord.gt2<-Wcoord.gt3<-Wcoord.gt4<-rep(0, nrow(data))

ll<-na.omit(data$Level)
ww<-na.omit(data$Wcoord)

Level.gt1[attributes(ll)[[1]]]<-Level.gt2[attributes(ll)[[1]]]<-Level.gt3[attributes(ll)[[1]]]<-Level.gt4[attributes(ll)[[1]]]<-Wcoord.gt1[attributes(ww)[[1]]]<-Wcoord.gt2[attributes(ww)[[1]]]<-Wcoord.gt3[attributes(ww)[[1]]]<-Wcoord.gt4[attributes(ww)[[1]]]<-NA

for (i in 1:nrow(data)){
if(data$Level[i]>1 & is.na(data$Level[i])==F) Level.gt1[i]<-1
  if(data$Level[i]>2 & is.na(data$Level[i])==F) Level.gt2[i]<-1
  if(data$Level[i]>3 & is.na(data$Level[i])==F) Level.gt3[i]<-1
if(data$Level[i]>4 & is.na(data$Level[i])==F) Level.gt4[i]<-1

if(data$Wcoord[i]>1 & is.na(data$Wcoord[i])==F) Wcoord.gt1[i]<-1
  if(data$Wcoord[i]>2 & is.na(data$Wcoord[i])==F) Wcoord.gt2[i]<-1
  if(data$Wcoord[i]>3 & is.na(data$Wcoord[i])==F) Wcoord.gt3[i]<-1
if(data$Wcoord[i]>4 & is.na(data$Wcoord[i])==F) Wcoord.gt4[i]<-1

}

data<-cbind(data,Level.gt1, Level.gt2, Level.gt3, Level.gt4, Wcoord.gt1, Wcoord.gt2, Wcoord.gt3, Wcoord.gt4)
data<-subset(data, year<2004) # no GWL data past 03

# interpolating intervening years for apphrfr1, noaff.largest.confed, pct.aff.largest.confed 

years<-seq(1948,2003)

form.a<- noaff.largest.confed ~ wbcode
form.b<- apphrfr1 ~ wbcode
form.c<- pct.aff.largest.confed ~wbcode

us.noaff<-cbind(unstack(data, form.a),years)
us.apphrfr1<-cbind(unstack(data,form.b),years)
us.pct.aff.largest.confed<-cbind(unstack(data,form.c),years)

missings<-c("ESP","NZL","FRA","PRT","IRL","GRC") # countries with no data whatever

us.noaff.interp<-us.noaff[,-22]
us.apphrfr1.interp<-us.apphrfr1[,-22]
us.pct.aff.largest.confed.interp<-us.pct.aff.largest.confed[,-22]


for(i in 1:ncol(us.apphrfr1)){
    if(colnames(us.noaff)[i] %in% missings)
    us.noaff.interp[,i]<-NA
    else
    us.noaff.interp[,i]<-na.approx(us.noaff[,i],along=years,na.rm=F)

    if(colnames(us.apphrfr1)[i] %in% missings)
    us.apphrfr1.interp[,i]<-NA
    else
    us.apphrfr1.interp[,i]<-na.approx(us.apphrfr1[,i],along=years,na.rm=F)

    if(colnames(us.pct.aff.largest.confed)[i] %in% missings)
    us.pct.aff.largest.confed.interp[,i]<-NA
    else
    us.pct.aff.largest.confed.interp[,i]<-na.approx(us.pct.aff.largest.confed[,i],along=years,na.rm=F)
    }

    
noaff.largest.confed.interp<-stack(us.noaff.interp[,-22])[,1]
apphrfr1.interp<-stack(us.apphrfr1.interp[,-22])[,1]
pct.aff.largest.confed.interp<-stack(us.pct.aff.largest.confed.interp[,-22])[,1]

data<-cbind(data,noaff.largest.confed.interp,apphrfr1.interp,pct.aff.largest.confed.interp)


# calculating eta^2 = log(H*a)

eta.sq <- log(noaff.largest.confed.interp * apphrfr1.interp)


# data matrix for analysis
data<-cbind(data, eta.sq)
DM<-na.omit(subset(data, select=c(id,wbcode,year,yr.trend, counter,trade,pop,enpp,left.gs,net.den, emp.contract.veto, eta.sq,no.confed,c.strike.fund.1, c.frac,federal, apphrfr1.interp, noaff.largest.confed.interp, pct.aff.largest.confed.interp)))


## TABLE 1 models

model.1<-geeglm(c.strike.fund.1 ~  eta.sq + pop + trade + net.den + no.confed
                + left.gs + enpp + as.factor(emp.contract.veto)+ yr.trend,
                id=wbcode, family=binomial(link=logit), data=DM)
model.1.test<-update(model.1, .~. - eta.sq)

model.2<-geeglm(c.strike.fund.1 ~  eta.sq + pop + trade + net.den + no.confed
                + left.gs + enpp + c.frac +as.factor(federal) +
                as.factor(emp.contract.veto), id=wbcode,
                family=binomial(link=logit), data=DM)
model.2.test<-update(model.2, .~. - eta.sq)


model.3<-geeglm(c.strike.fund.1 ~  eta.sq + pop + trade + net.den + no.confed
                + left.gs + enpp + c.frac +as.factor(federal)
                + as.factor(emp.contract.veto) + yr.trend,
                id=wbcode, family=binomial(link=logit), data=DM)
model.3.test<-update(model.3, .~. - eta.sq)


pred.strk.fnd<-fitted(model.2) # needed for predicted strike fund.
DM.plus<-cbind(DM, pred.strk.fnd) #  use below

# WSEV calcs for table 1 models
X.1<-model.matrix(model.1)
yrs<-data$year[as.numeric(row.names(X.1))]
m<-10

X.2<-model.matrix(model.2)
X.3<-model.matrix(model.3)


strk.fnd.1.var<-window.var(y=model.1$y, X=X.1, fit = model.1,
                             tvec=yrs, m=m)
se.1<-strk.fnd.1.var$window.1

strk.fnd.2.var<-window.var(y=model.2$y, X=X.2, fit = model.2,
  tvec=yrs, m=m)
se.2<-strk.fnd.2.var$window.1

strk.fnd.3.var<-window.var(y=model.3$y, X=X.3, fit = model.3,
  tvec=yrs, m=m)
se.3<-strk.fnd.3.var$window.1


model.1[[29]]<-se.1
names(model.1)[[29]]<-"se"

model.2[[29]]<-se.2
names(model.2)[[29]]<-"se"

model.3[[29]]<-se.3
names(model.3)[[29]]<-"se"

table.edit(model.1,model.2, model.3, se="robust", stars=1,lev=.0001,  order="rl") #hacked version of apsrtable


## robustness fact mentioned in fn 13
d2<-na.omit(subset(data, select=c(id,wbcode,year,yr.trend, counter,trade,pop,enpp,left.gs,net.den, emp.contract.veto, eta.sq,no.confed,Cfauthority, c.frac,federal, apphrfr1.interp, noaff.largest.confed.interp, pct.aff.largest.confed.interp)))
summary(geeglm(Cfauthority ~  eta.sq + pop + trade + net.den + no.confed
                + left.gs + enpp + as.factor(emp.contract.veto)+ yr.trend, id=wbcode,data=d2))

## separation plots (figure 1)
source("separationplot.R",local=T) #contact Brian Greenhill for separation plot code


main.1<-"model 1 excluding"~ eta^2
separationplot(pred=as.vector(model.1$fitted), actual=model.1$y, heading = "model 1", file="sp1.pdf")
separationplot(pre=as.vector(model.1.test$fitted), actual=model.1$y, heading = main.1, file="sp1test.pdf")

main.2<-"model 2 excluding"~ eta^2
sp2<-separationplot(pred=as.vector(model.2$fitted), actual=model.2$y, heading = "model 2", file="sp2.pdf")
sp2.test<-separationplot(pre=as.vector(model.2.test$fitted), actual=model.2$y, heading = main.2, file="sp2test.pdf")

main.3<-"model 3 excluding"~ eta^2
sp3<-separationplot(pred=as.vector(model.3$fitted), actual=model.3$y, heading = "model 3", file="sp3.pdf")
sp3.test<-separationplot(pre=as.vector(model.3.test$fitted), actual=model.3$y, heading = main.3, file="sp3test.pdf")


## Ladder plot (fig 2) using model 3

pred.means<-apply(X.3,2,mean)                                     
pred.means[6]<-2 #integer value for no.confed
pred.means[c(10,11)]<-0 # modal values for emp. veto & federalism
pred.means<-pred.means[-1]  # get rid of constant

pred.iq<-apply(X.3,2,quantile,c(.25,.75)) #interquartile range
pred.iq<-pred.iq[,-1]

b.sim<-mvrnorm(1000,model.3$coef,se.3) #drawing coefficient vectors from posterior

lamda.1q<-lamda.3q<-matrix(0,nrow=length(pred.means),ncol=1000)
for(i in 1:length(pred.means)){
x.1q<-x.3q<-pred.means
x.1q[i]<-pred.iq[1,i]
x.3q[i]<-pred.iq[2,i]
x.1q<-c(1,x.1q)
x.3q<-c(1,x.3q)
  for(j in 1:1000){
    lamda.1q[i,j]<-exp(x.1q%*%b.sim[j,])/(1+exp(x.1q%*%b.sim[j,]))
    lamda.3q[i,j]<-exp(x.3q%*%b.sim[j,])/(1+exp(x.3q%*%b.sim[j,]))
    }
 }



diff<-lamda.3q-lamda.1q
diff.mean<-apply(diff,1,mean)
diff.ci<-apply(diff,1,quantile,c(0.025,.975))
diff.mean<-diff.mean[-length(diff.mean)] # get rid of year trend
diff.ci<-diff.ci[,-(length(diff.mean)+1)] # get rid of year trend

ivs.logit<-c(expression(eta^2),"population", "trade", "density", "no. confed", "left gov't","ENPP", "cultural frac.",
  "federal","emp. cntrct veto")

pdf("strkfnd.pdf")
par(mar=c(2,6,4,2))
plot(0,0, xlim=c(-1,1), ylim = c(10,110), type="n",xlab="",ylab="", yaxt="n",main="Changes in predicted probabilities", bty="n")
abline(v=0, lty=2, col=grey(.5))
mtext(ivs.logit, side=2, at=c(105,95,85,75,65,55,45,35,25,15),font=1, padj=0, las=2)
points(diff.mean,c(105,95,85,75,65,55,45,35,25,15))
segments(x0=diff.ci[1,], y0=c(105,95,85,75,65,55,45,35,25,15), x1=diff.ci[2,], y1=c(105,95,85,75,65,55,45,35,25,15), lwd=2)
dev.off()

#### Bargaining coordination models

# getting predicted prob of strike fund
p.strike.fund<-p.emp.veto<-rep(NA, nrow(data))
p.strike.fund[match(DM.plus$id, data$id)]<-DM.plus$pred.strk.fnd
data<-cbind(data,p.strike.fund)

#setting up data matrix for analysis
DM.Level<-na.omit(subset(data, select=c(id,wbcode,year,yr.trend, counter,trade,pop,enpp,left.gs,net.den,
                                    emp.contract.veto, eta.sq,no.confed,c.strike.fund.1, c.frac,federal,
                                    Level.gt1, Level.gt2, Level.gt3, Level.gt4,
                                    noaff.largest.confed.interp, pct.aff.largest.confed.interp,
                                    p.strike.fund)))

DM.Wcoord<-na.omit(subset(data, select=c(id,wbcode,year,yr.trend, counter,trade,pop,enpp,left.gs,net.den,
                                    emp.contract.veto, eta.sq,no.confed,c.strike.fund.1, c.frac,federal,
                                    Wcoord.gt1, Wcoord.gt2, Wcoord.gt3, Wcoord.gt4,
                                    noaff.largest.confed.interp, pct.aff.largest.confed.interp,
                                    p.strike.fund)))

# Visser Wcoord variable
wcoord.gt1.fit<-geeglm(Wcoord.gt1~ p.strike.fund + eta.sq + pop + trade + net.den + no.confed + left.gs
                   + federal + yr.trend ,
                  data=DM.Wcoord, id=wbcode, family=binomial) # no c.frac, enpp, or emp.contract.veto b/c perfect separation. 
wcoord.gt2.fit<-geeglm(Wcoord.gt2~ p.strike.fund + eta.sq +pop + trade + net.den + no.confed + left.gs
                  + c.frac + federal +  yr.trend ,
                  data=DM.Wcoord, id=wbcode, family=binomial) # no enpp, emp.contract.veto b/c sigular WSEV 
wcoord.gt3.fit<-geeglm(Wcoord.gt3~ p.strike.fund + eta.sq +pop + trade + net.den + no.confed + left.gs
                  + enpp   + c.frac + federal + emp.contract.veto + yr.trend ,
                  data=DM.Wcoord, id=wbcode, family=binomial)
wcoord.gt4.fit<-geeglm(Wcoord.gt4~ p.strike.fund + eta.sq + pop + trade + net.den + no.confed + left.gs
                  + enpp   + c.frac + federal + emp.contract.veto + yr.trend ,
                  data=DM.Wcoord, id=wbcode, family=binomial)

## calc WSEV
                  
X.gt1<-model.matrix(wcoord.gt1.fit)
X.gt2<-model.matrix(wcoord.gt2.fit)
X.gt3<-model.matrix(wcoord.gt3.fit)
X.gt4<-model.matrix(wcoord.gt4.fit)

yrs<-data$year[as.numeric(row.names(X.gt1))]
m<-10

gt1.var<-window.var(y=wcoord.gt1.fit$y, X=X.gt1, fit = wcoord.gt1.fit,
  tvec=yrs, m=m)
se.gt1<-gt1.var$window.1

gt2.var<-window.var(y=wcoord.gt2.fit$y, X=X.gt2, fit = wcoord.gt2.fit,
  tvec=yrs, m=m)
se.gt2<-gt2.var$window.1

gt3.var<-window.var(y=wcoord.gt3.fit$y, X=X.gt3, fit = wcoord.gt3.fit,
  tvec=yrs, m=m)
se.gt3<-gt3.var$window.1

gt4.var<-window.var(y=wcoord.gt4.fit$y, X=X.gt4, fit = wcoord.gt4.fit,
  tvec=yrs, m=m)
se.gt4<-gt4.var$window.1


wcoord.gt1.fit[[29]]<-se.gt1
names(wcoord.gt1.fit)[[29]]<-"se"

wcoord.gt2.fit[[29]]<-se.gt2
names(wcoord.gt2.fit)[[29]]<-"se"

wcoord.gt3.fit[[29]]<-se.gt3
names(wcoord.gt3.fit)[[29]]<-"se"

wcoord.gt4.fit[[29]]<-se.gt4
names(wcoord.gt4.fit)[[29]]<-"se"



table.edit(wcoord.gt1.fit, wcoord.gt2.fit, wcoord.gt3.fit, wcoord.gt4.fit, se="robust", stars=1,lev=.00001,  order="rl")

## Ladder plot (fig 3)

X1<-X.gt1
X2<-X.gt2
X3<-X.gt3 #X.gt.3 = X.gt4

#X1 is missing enpp, c.frac, and emp.contract.veto
#X2 is missing enpp and emp.contract.veto

pred.means<-apply(X1,2,mean)
pred.means<-pred.means[-1]

pred.means2<-apply(X2,2,mean)
pred.means2<-pred.means2[-1]

pred.means3<-apply(X3,2,mean)
pred.means3<-pred.means3[-1]

pred.means[6]<-pred.means2[6]<-pred.means3[6]<-2  # no. confed
pred.means[8]<-0
pred.means2[9]<-0
pred.means3[10]<-0  # fed


pred.iq<-apply(X1,2,quantile,c(.25,.75))
pred.iq<-pred.iq[,-1]

pred.iq2<-apply(X2,2,quantile,c(.25,.75))
pred.iq2<-pred.iq2[,-1]

pred.iq3<-apply(X3,2,quantile,c(.25,.75))
pred.iq3<-pred.iq3[,-1]

b1.sim<-mvrnorm(1000,wcoord.gt1.fit$coef,se.gt1)
b2.sim<-mvrnorm(1000,wcoord.gt2.fit$coef,se.gt2)
b3.sim<-mvrnorm(1000,wcoord.gt3.fit$coef,se.gt3)
b4.sim<-mvrnorm(1000,wcoord.gt4.fit$coef,se.gt4)

lamda1.1q<-lamda1.3q<-matrix(0,nrow=length(pred.means),ncol=1000)
lamda2.1q<-lamda2.3q<-matrix(0,nrow=length(pred.means2),ncol=1000)
lamda4.1q<-lamda4.3q<-lamda3.1q<-lamda3.3q<-matrix(0,nrow=length(pred.means3),
                                                   ncol=1000)


for(i in 1:length(pred.means)){
x.1q<-x.3q<-pred.means
x.1q[i]<-pred.iq[1,i]
x.3q[i]<-pred.iq[2,i]
x.1q<-c(1,x.1q)
x.3q<-c(1,x.3q)
  for(j in 1:1000){
    lamda1.1q[i,j]<-exp(x.1q%*%b1.sim[j,])/(1+exp(x.1q%*%b1.sim[j,]))
    lamda1.3q[i,j]<-exp(x.3q%*%b1.sim[j,])/(1+exp(x.3q%*%b1.sim[j,]))
  }}

for(i in 1:length(pred.means2)){
x.1q<-x.3q<-pred.means2
x.1q[i]<-pred.iq2[1,i]
x.3q[i]<-pred.iq2[2,i]
x.1q<-c(1,x.1q)
x.3q<-c(1,x.3q)
  for(j in 1:1000){
    lamda2.1q[i,j]<-exp(x.1q%*%b2.sim[j,])/(1+exp(x.1q%*%b2.sim[j,]))
    lamda2.3q[i,j]<-exp(x.3q%*%b2.sim[j,])/(1+exp(x.3q%*%b2.sim[j,]))
    }
 }

for(i in 1:length(pred.means3)){
x.1q<-x.3q<-pred.means3
x.1q[i]<-pred.iq3[1,i]
x.3q[i]<-pred.iq3[2,i]
x.1q<-c(1,x.1q)
x.3q<-c(1,x.3q)
  for(j in 1:1000){
    lamda3.1q[i,j]<-exp(x.1q%*%b3.sim[j,])/(1+exp(x.1q%*%b3.sim[j,]))
    lamda3.3q[i,j]<-exp(x.3q%*%b3.sim[j,])/(1+exp(x.3q%*%b3.sim[j,]))

    lamda4.1q[i,j]<-exp(x.1q%*%b4.sim[j,])/(1+exp(x.1q%*%b4.sim[j,]))
    lamda4.3q[i,j]<-exp(x.3q%*%b4.sim[j,])/(1+exp(x.3q%*%b4.sim[j,]))
  }
}
diff1<-lamda1.3q-lamda1.1q
diff2<-lamda2.3q-lamda2.1q
diff3<-lamda3.3q-lamda3.1q
diff4<-lamda4.3q-lamda4.1q


diff1.mean<-apply(diff1,1,mean)
diff2.mean<-apply(diff2,1,mean)
diff3.mean<-apply(diff3,1,mean)
diff4.mean<-apply(diff4,1,mean)

diff1.q<-apply(diff1,1,quantile,c(.975,.025))
diff2.q<-apply(diff2,1,quantile,c(.975,.025))
diff3.q<-apply(diff3,1,quantile,c(.975,.025))
diff4.q<-apply(diff4,1,quantile,c(.975,.025))

diff1.q<-diff1.q[,-length(diff1.mean)]
diff2.q<-diff2.q[,-length(diff2.mean)]
diff3.q<-diff3.q[,-length(diff3.mean)]
diff4.q<-diff4.q[,-length(diff4.mean)]

diff1.mean<-diff1.mean[-length(diff1.mean)]
diff2.mean<-diff2.mean[-length(diff2.mean)]
diff3.mean<-diff3.mean[-length(diff3.mean)]
diff4.mean<-diff4.mean[-length(diff4.mean)]

diff1.m.plot<-rep(2,length(diff4.mean))
diff1.q.plot<-matrix(0,nrow=2, ncol=length(diff4.mean))
diff1.m.plot[1:7]<-diff1.mean[1:7]
diff1.m.plot[10]<-diff1.mean[8]
diff1.q.plot[,1:7]<-diff1.q[,1:7]
diff1.q.plot[,10]<-diff1.q[,8]

diff2.m.plot<-rep(2,length(diff4.mean))
diff2.q.plot<-matrix(0,nrow=2, ncol=length(diff4.mean))
diff2.m.plot[1:7]<-diff2.mean[1:7]
diff2.m.plot[9:10]<-diff2.mean[8:9]
diff2.q.plot[,1:7]<-diff2.q[,1:7]
diff2.q.plot[,9:10]<-diff2.q[,8:9]

ivs.slogit<-c("Strike fund",expression(eta^2),"Population", "Trade","Density",
              "No. confed", "Left gov't", "ENPP", "Cultural frac", "Federal",
              "Emp cntrct veto")

pdf("coord.pdf")
par(mar=c(2,6,4,2))
plot(0,0, xlim=c(-1,1), ylim = c(3,132), type="n",xlab="",ylab="", yaxt="n",main="Changes in predicted probabilities", bty="n")
abline(v=0, lty=2, col=grey(.5))
mtext(ivs.slogit, side=2, at=c(132-6,119-6,107-6,95-6,83-6,71-6,59-6,47-6,35-6,23-6,11-6),font=1, padj=0, las=2)

segments(x0=diff1.q.plot[1,], y0=c(131.5,118.5,106.5,94.5,82.5,70.5,58.5,46.5,34.5,22.5,10.5),
         x1=diff1.q.plot[2,], y1=c(131.5,118.5,106.5,94.5,82.5,70.5,58.5,46.5,34.5,22.5,10.5),
         lwd=2, col=grey(0.5))

segments(x0=diff2.q.plot[1,], y0=c(129,116,104,92,80,68,56,44,32,20,8),
         x1=diff2.q.plot[2,], y1=c(129,116,104,92,80,68,56,44,32,20,8),
         lwd=2, col=grey(0.5))

segments(x0=diff3.q[1,], y0=c(126.5,113.5,101.5,89.5,77.5,65.5,53.5,41.5,29.5,17.5,5.5),
         x1=diff3.q[2,], y1=c(126.5,113.5,101.5,89.5,77.5,65.5,53.5,41.5,29.5,17.5,5.5),
         lwd=2, col=grey(0.5))

segments(x0=diff4.q[1,], y0=c(124,111,99,87,75,63,51,39,27,15,3),
         x1=diff4.q[2,], y1=c(124,111,99,87,75,63,51,39,27,15,3),
         lwd=2, col=grey(0.5))


points(diff1.m.plot,c(131.5,118.5,106.5,94.5,82.5,70.5,58.5,46.5,34.5,22.5,10.5), pch=1)
points(diff2.m.plot,c(129,116,104,92,80,68,56,44,32,20,8), pch=17)
points(diff3.mean,c(126.5,113.5,101.5,89.5,77.5,65.5,53.5,41.5,29.5,17.5,5.5), pch=22)
points(diff4.mean,c(124,111,99,87,75,63,51,39,27,15,3), pch=19)

abline(h=c(120,108,96,84,72,60,48,36,24,12), col=grey(.3), lwd=.5)

dev.off()

## Robustness check for Level variable mentioned in fn 23
#Visser Level variable
level.gt1.fit<-geeglm(Level.gt1~ p.strike.fund + eta.sq + pop + trade + net.den + no.confed + left.gs
                   + federal + yr.trend ,
                  data=DM.Level, id=wbcode, family=binomial) # no emp.contract.veto b/c singular WSEV
level.gt2.fit<-geeglm(Level.gt2~ p.strike.fund + eta.sq + pop + trade + net.den + no.confed + left.gs
                  + enpp  + federal + emp.contract.veto +  yr.trend ,
                  data=DM.Level, id=wbcode, family=binomial) 
level.gt3.fit<-geeglm(Level.gt3~ p.strike.fund + eta.sq + pop + trade + net.den + no.confed + left.gs
                  + enpp   + c.frac + federal + emp.contract.veto + yr.trend ,
                  data=DM.Level, id=wbcode, family=binomial)
level.gt4.fit<-geeglm(Level.gt4~ p.strike.fund + eta.sq + pop + trade + net.den + no.confed + left.gs
                  + enpp + emp.contract.veto + yr.trend ,
                  data=DM.Level, id=wbcode, family=binomial) # no federal, c.frac b/c WSEV singular


#WSEV for level

X.gt1<-model.matrix(level.gt1.fit)
X.gt2<-model.matrix(level.gt2.fit)
X.gt3<-model.matrix(level.gt3.fit)
X.gt4<-model.matrix(level.gt4.fit)

yrs<-data$year[as.numeric(row.names(X.gt1))]
m<-10

gt1.var<-window.var(y=level.gt1.fit$y, X=X.gt1, fit = level.gt1.fit,
  tvec=yrs, m=m)
se.gt1<-gt1.var$window.1

gt2.var<-window.var(y=level.gt2.fit$y, X=X.gt2, fit = level.gt2.fit,
  tvec=yrs, m=m)
se.gt2<-gt2.var$window.1

gt3.var<-window.var(y=level.gt3.fit$y, X=X.gt3, fit = level.gt3.fit,
  tvec=yrs, m=m)
se.gt3<-gt3.var$window.1

gt4.var<-window.var(y=level.gt4.fit$y, X=X.gt4, fit = level.gt4.fit,
  tvec=yrs, m=m)
se.gt4<-gt4.var$window.1


level.gt1.fit[[29]]<-se.gt1
names(level.gt1.fit)[[29]]<-"se"

level.gt2.fit[[29]]<-se.gt2
names(level.gt2.fit)[[29]]<-"se"

level.gt3.fit[[29]]<-se.gt3
names(level.gt3.fit)[[29]]<-"se"

level.gt4.fit[[29]]<-se.gt4
names(level.gt4.fit)[[29]]<-"se"

table.edit(level.gt1.fit, level.gt2.fit, level.gt3.fit, level.gt4.fit, se="robust", stars=1,lev=.00001,  order="rl")

### END analysis code
#### END
